home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / SCIENTIF / H381.ZIP / GSRC208A.ZIP / NEWTON.C < prev    next >
C/C++ Source or Header  |  1993-07-16  |  4KB  |  109 lines
  1. #include "copyleft.h"
  2.  
  3. /*
  4.     GEPASI - a simulator of metabolic pathways and other dynamical systems
  5.     Copyright (C) 1989, 1992  Pedro Mendes
  6. */
  7.  
  8. /*************************************/
  9. /*                                   */
  10. /*      steady-state solution by     */
  11. /*      the damped Newton method     */
  12. /*                                   */
  13. /*        Zortech C/C++ 3.0 r4       */
  14. /*          MICROSOFT C 6.00         */
  15. /*          Visual C/C++ 1.0         */
  16. /*           QuickC/WIN 1.0          */
  17. /*             ULTRIX cc             */
  18. /*              GNU gcc              */
  19. /*                                   */
  20. /*   (include here compilers that    */
  21. /*   compiled GEPASI successfully)   */
  22. /*                                   */
  23. /*************************************/
  24.  
  25.  
  26. #include <stdio.h>
  27. #include <math.h>
  28. #include <string.h>
  29. #include <stdlib.h>
  30. #include "globals.h"                       /* global symbols                 */
  31. #include "globvar.h"                       /* extern global variable         */
  32. #include "matrix.h"                        /* matrix manipulation            */
  33. #include "rates.h"                         /* rates and jacobian calculation */
  34. #include "datab.h"                           /* database structures             */
  35.  
  36. #define N_OK 0
  37. #define N_LMIN 1
  38. #define N_NOCONV 2
  39. #define N_JSING 3
  40.  
  41. int newton( void )
  42. {
  43.  register int m, i, j, k;
  44.  unsigned long df;                         /* dumping factor (power of two)  */
  45.  static double xpss[MAX_MET],              /* x at iteration m               */
  46.                ma[MAX_MET][MAX_MET],       /* inverse of jacobian matrix     */
  47.                mb[MAX_MET][MAX_MET],       /* copy of jacobian matrix        */
  48.                h[MAX_MET], h2[MAX_MET],    /* increment and damped increment */
  49.                re,nre;                     /* residual error & new res.err.  */
  50.  
  51.  
  52.  for(j=0;j<totmet;j++)                     /* make x[i] the first guess      */
  53.   xpss[j] = xss[j] = x[j];
  54.                                            /* initial values of f(x) and re  */
  55.  calcrates( xpss );                        /* f(x)                           */
  56.  for(j=0, re = (double) 0; j<nmetab;j++)   /* ||f(x)||2 residual error       */
  57.   re += rate[j] * rate[j];
  58.  re = sqrt( re );
  59.  if ( re < options.hrcz ) return N_OK;     /* if already in s.s. job is done!*/
  60.  for(m=0;m<newtlim;m++)                    /* iterate "until satisfied"      */
  61.  {
  62.   if (options.debug) printf("\nnewton() - iteration %2d, ", m);
  63.   calcjacob( xpss );                       /* f'(x) jacobian                 */
  64.   memcpy( mb, jacob, MAX_MET * MAX_MET * sizeof(double) );
  65.   lu_inverse( (double(*)[MAX_MET][MAX_MET])mb,
  66.               (double(*)[MAX_MET][MAX_MET])ma,
  67.               indmet );                    /* inv( f'(x) )                   */
  68.   multmtxv( ma, rate, h,
  69.             indmet, indmet, indmet );      /* h = inv(f'(x)) . f(x)          */
  70.   nre = 2 * re;                            /* just to start with nre > re    */
  71.   for( i=0, df=1; (i<32) && (nre>re) ; i++ )
  72.   {
  73.    if (options.debug) printf(".");
  74.    if (i) df *= 2;                         /* 2 raised to the i power        */
  75.    for( j=0; j<indmet; j++ )
  76.     h2[j] = h[j]/df;
  77.    addvctsm( xpss, h2, -1.0, xss, indmet );/* xss = xpss - h2                */
  78.    /* update the conserved moieties            */
  79.    for( j=indmet; j<nmetab; j++)
  80.    {
  81.     xss[j] = moiety[j];
  82.     for( k=0; k<indmet; k++)
  83.      xss[j] -= ld[j][k] * xss[k];
  84.    }
  85.    calcrates( xss );                        /* f(x)                           */
  86.    for(j=0, nre=0.0; j<indmet; j++)            /* ||f(x)||2 residual error       */
  87.     nre += rate[j] * rate[j];
  88.    nre = sqrt( nre );
  89.   }
  90.   re = nre;                                    /* keep nre for next step         */
  91.   if (i==32)
  92.   {
  93.    if ( re < options.hrcz )  return N_OK;   /* this is a zero, after all!     */
  94.    else                                        /* stuck in a local minima          */
  95.    {
  96.     if (options.debug) printf(" local minima ");
  97.     return N_LMIN;
  98.    }
  99.   }
  100.   if (options.debug) printf(" re=% .4le", re );
  101.   for(j=0; j<nmetab; j++)
  102.    xpss[j] = xss[j];
  103.   if ( re < options.hrcz )
  104.    return N_OK;
  105.  }
  106.  if (options.debug) printf(" no solution ");
  107.  return N_NOCONV;                          /* no convergence                  */
  108. }
  109.